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Abstract 

New numerical method to calculate thermodynmic Bethe ansatz equations is proposed based 
on Newton's method. Thermodynamic quantities of one-dimensional Hubbard model is numeri- 
cally calculated and compared with high temperature expansion and numerical results of quantum 
transfer matrix method by Jiittner, Kliimper and Suzuki. The coincidence is surprisingly good. 
We get high-temperature expansion of grand potential up to (3^. 

PACS numbers: 71. 27. -Fa, 05.30.-d, 05.30.Fk 
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I. INTRODUCTION 



Many years ago, one of the authors (MT) proposed thermodynamic Bethe ansatz (TBA) 
equations for one-dimensional Hubbard model [2]. In this theory several kind of strings are 
assumed and it is widely believed that this set of equations give the exact thermodynamic 
quantities of this model. Low-temperature thermodynamics were investigated by MT [3] 
and actual numerical calculations at finite temperature were done by Kawakami, Usuki and 
Okiji [4]. Essler, Korepin and Schoutens [5] counted the number of states by the single 
k excitations, A strings and k — A strings. They found that total number of these Bethe 
ansatz states and their relatives is 4^", where Na is the length of the systems. This implies 
that the Bethe ansatz can give all eienstates and eigenvalues. However some physicists are 
still skeptical for this theory [8]. Recently Charret et al [7] did the numerical calculation 
of this equation and concluded that it does not coincide with high temperature expansion 
and quantum transfer matrix(QTM) method by Jiittner, Kliimper and Suzuki [6]. In this 
paper we give a practical method to calculate numerically TBA equations which has infinite 
unknown functions. Numerical results completely coincide with those of QTM and HTE. 
Charret et al's numerical calculation of TBA equations is wrong. The Hubbard Hamiltonian 
is 



Here cj^ and Cj^j are creation and annihilation operators of an electron at site j. < ij > means 
that sites i and j are nearest neighbors. Na is the number of atoms. We put t > 0, [/ > 0. 
Thermodynamic potential per site g at temperature T is determined by 



n{tjj,A,h) 






g = cq — A — T 
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Here eo,/3o(^) and cro(A) are energy per cite, distribution functions of A;'s and A's at T 
h = U /2 — A = 0^ (half-filled, zero-field ground state). 



oo 







(jj 



[1 + exp{2U'uj))' 



ao(A)= [ siA-smk)^, (4) 

f 

po{k) = h COS A; / ai(A — sin A;)(To(A)c/A, (5) 



and 



C{k) and ?7i(A) are hole-particle ratios of k excitations and single A excitations. These are 
determined by thermodynamic Bethe ansatz equations for k excitations, A strings and k — A 
strings; 

Ko{k) 



In C{k) 



T 



+ £.A,„A-s„U.)ln(i±|||). (7) 

ln?7i(A) = s * ln(f + r]2{A)) 

/TV 
s{A - smk)ln{l + C-\k)) cos kdk, (8) 
-TT 

ln?7;(A) = s*ln(f + ?7^(A)) 

/TT 
s(A - sinA;)ln(f + C(A;))cosA;(/A;, (9) 

lnr7,(A) = s * ln{(f + r],-r{A)){l + r7,+i(A))}, 

J > 2, (10) 
lnry;(A) = s * ln{(l + ry;_,(A))(l + /^^^^(A))}, 

J > 2, (11) 



ln?7„(A) 2h , , 

lim = — , (12) 

ra->oo n i 

li„. = ^i^. (13) 

ra->oo n r 



Here s * /(A) = s(A - A')f{A')dA' and /€o(A;) is defined by 

Ko[k) = — 2t cos k 



At / c/As(A - sinA;)3f?Vl - (A - U'lf. (14) 



OO 



Details of derivation are given in references [1, 2, 9, 10]. Noting that ln(l + C"^) = ln(l + 
C) — InC in (8) and substituting (7) we have 



In Ref. [2] this set of equations were solved analytically in the limits T — > 0, t — > 0, and 
[/ — > and coincided with known exact results. In a recent paper Charret et al [7] solved 
numerically this set of equations at high temperature and argued that there is discrepancy 
from high temperature expansion and numerical results of Jiittner, Kliimper and Suzuki 
equations [6]. We recalculate the same quantities in this region and find that the results 
coincide with high temperature expansion and JKS equations in high accuracy. In III we 
review t expansion for one-dimensional Hubbard model by the conventional linked cluster 
expansion. From the t expansion, we can derive (3 expansion of —g(3 for the ID Hubbard 
model up to (3^. Expansions of susceptibility and specific heat are obtained. In Appendix 
A, we can perform t expansion of TEA equations. The results coincide with the cluster 
expansion up to the second order. We expect that the higher terms also coincide. 



II. TRUNCATION OF TBA EQUATIONS TO FINITE UNKNOWN FUNCTIONS 



ln?7i(A) 





(15) 



As an approximation we replace s(A) by |<^(A) at j > ric in equations 
Then we get the difference equations 




r7,(Af = (l + r7,_i(A))(l + r7,+i(A)), 




(16) 



This approximation is reasonable because functions ?7j(A) and ?7j(A) vary very slowly at 
sufficiently large j and s(A)c/A = 1/2. General solutions of these difference equations 



are 
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From the conditions (12) and (13), parameters a and b must be h/T and [U/2 — A)/T. Then 
1 + ?7„^+i(A) and 1 + ri'^^_^_^[A) are represented by 

(cosh ^y/l + r]^^{A) + ^l + sinh2^(l + r7„,(A)))', 

(cosh iyrr^^A) + ^l + sinh2^(l + <^(A)))', 

h' = U/2 - A. (18) 

Thus integral equations with infinite unknown functions are approximated by those with 
2nc + 1 unknowns In ?7^^(A), ln?7^(A), In C{k)^ In ?7i(A), ln?7„^(A). Then equations to be 
solved are 



Zi = s * ln[(l + exp 2;2)(cosh U2 y 1 + exp Zi -\- y 1 -\- sinh ^2(1 + exp Zi)) 



Z4 



J - ln(l + exp2;j_i)(l + exp^r^+i), j = 2, ...,12^ - 1, 



Z.n 



^ = s * ln(l + exp 2;„^_i) — / s(A — sin A;) ln(l + exp 2;„^_l_i) cos A;c/A;, 



^+1 = ui^o + / c/As(A — sin A;) ln(- ~'~ P 



-00 



1 + exp Zn, + 2 

/TT 
s(A — sin k) ln(l + exp Zn^+i) cos kdk, 
■TT 

= s * ln(l + exp2;j_i)(l + exp2;j+i), j = ric + 3, ...,2^^, 



-22ne+i = * ln[(l + exp 2;2„ J( cosh usy/l + exp2;2„,+i + y 1 + sinh 1/3(1 + exp2;2„,+i)) ]. 

(19) 

We introduce three thermodynamic parameters: 

ui = l/r, U2 = {U/2 - A)/T, U3 = h/T. 
For actual numerical calculations we choose L discrete points of k and A as follows: 



kj = T\:{j — 1/2)/ L, Aj = s'mqjill + 



COS^ 



q, = 7T{j-l/2)/{2L), j = l,...,I. (20) 

Here function A = sin qWl + {W / cos qY is the inverse function of P ^{1 -{t- U'lfY^l^dt. 
We think that this change of parameters is reasonable because the change of functions is very 
slow at large A. For very big U' this function behaves as U' tan q and for small U' it behaves 



as sin q. Unknown functions are represented by vectors with length L and integration kernels 
are represented hy L X L matrices. 

Usually this kind of non-linear equations is calculated by successive iterations, which is 
called Kepler's method. In the solution of TBA equations we need to repeat several tens or 
several hundreds times of iterations to get a good convergence. 

So here we propose to use Newton's method. Consider a coupled non-linear equations: 

X,-F,(Xi,X2,...,Xa.) = 0, j = l,...,iV. (21) 
For approximate vectors xj'^ assume that we have deviations A^-: 

xf^-F,{x['\xi'\...,xii^) = A,. (22) 
In Kepler's method next approximation is 

xf+')=xf + A,. (23) 

In Newton's method we put 

Xf+^) = xf + e„ (24) 
where is the solution of linear equation 

Z^iki " ^ ^ 

J •' 

This method is much faster than Kepler's method. But we must solve linear equations with 
N X N matrix. In our TBA problem, N is [2nc + 1)L. This large matrix is block tridiagonal. 
Regarding L X L blocks as a number we can solve this set of linear equations. We need only 
5-6 times of iterations at most to get sufficient convergence |Aj| < 10~®. We can get the 
thermodynamic potential through equation (2): 

g , U 

po{k) ln(l + exp Zn,+i {k))dk 

ao{A) ln(l + exp z„,+2(A))(/A. (26) 
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To get the first order thermodynamic quantities like magnetization (m), electron density 
[n) and entropy we need to calculate d{g /T) / dui^ d{g /T) / du2^ d{g jT) j du^. 

d^j; = (eo - + 

po{k)— TTTTC!^ 

The equations for diZa is a linear equation which has the same homogeneous term with 
that in Newton's method. Inhomogeneous terms are calculated from Za- Therefore we can 
calculate these quantities by one operation of linear calculation, 

e = d,{g/T) + Ad2{g/T), 

n = d2{g/T), m = daig/T), 

entropy = ui[e — g) — u^m + (u2 — Uui/2)n. (28) 

To calculate the second order thermodynamic quantities such as specific heat, suscepti- 
bility and compressibility we need 3x3 tensor didj{g /T). As this is a symmetric tensor, we 
need to calculate six components. We consider equations for 

= d,d,z„ + {d,z„){d,z„), I < J. (29) 

1 + exp Za 



Tensor is given by 



7f = - Po{k)— -7--' 

I J 1 + exp(-2;„^+i(A;)) 



f^'^'^ (A) 

^o(A)- ■^";+'^ ' ,^,, dA. (30) 

V ^i + exp(-z„,+2(A)) ^ ^ 

The inhomogeneous terms are calculated from Za and diZa- So we can calculate thermody- 
namic quantities from 10 quantities g/T and its derivatives for given temperature, magnetic 
field and chemical potential after several times of linear equation solving. Specific heat c 
in constant n and m, susceptibility x in constant n and temperature, compressibility k in 
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constant temperature and m/n are given by 



c = uKDu + 
X = uADss - 



23 N 



K = Ui 



D 



D22 

(nDi3 + mD22)D23 



22 



nD33 + mD 
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A, = d,d,ig/T). 



(31) 



If we want to calculate thermodynamic quantities at fixed electron density, we can use 
Newton's method again. Using the compressibility we reach target density after few times 
of iterations. 



III. i-EXPANSION AND /3(= l/T) EXPANSION 

The t expansion of thermodynamic potential up to t"^ has been done for single band 
Hubbard model by Kubo [11] and Liu [12]. For one-dimensional model their expansion 
becomes 

-g(3 = iog[e] + m'^G, + m\^G2 + ^03 - + omf), (32) 

where 

e = l + 2cosh(/3/i)e^^ + e^(2A-c/)^ 



e 



G, = cosH(3h)e^\l + e^(2^-^)) + ^(1 

pU 

03 = -ee^^ cosh((3h)(l + e/^(2^-^)) + ^(1 + e/^(2^-^) - 2 cosh(/3/.)e^(^-^)) 
6 pU 

+ |^(2cosh(/3/.)e^^(l - 2e-n - (2 - e-^^)(l + e^(^^-^))) + ^(1 - e"^^). 

(33) 

Gi term is the second order term of t expansion. G2, G3 and Gf terms are the fourth order. It 
is expected that —g(3 is expanded by (3h^ (3h\ fiu\ where we put h' = U /2 — u' = U /4: 

-gf3 = J2 A^,,^,,^,,^,{(3tr{(3hr{m^^^ (34) 

ni ,n2 ,n3 ,n4 >0 



From equation (32) and (33) we can calculate all coefficients An^,^^,^^^,^^ at n\ < 6. On the 
other hand we have 



{jstf iptr m! 

16 144 



= ln4 + ^ - ^ + ^ - Om)% (35) 



when U = h = h' = 0. Then we have A6,o,o,o = 1/144 and ^7^0,0,0 = 0. In this way, we can 
determine An^,^^,^^^,^^ at ni + ^2 + ^3 + ^4 < 7 and obtain /3 expansion of the grand potential: 

g = + {h' - u') - ^{2t' + 2u" + h' + h") - ^{h' - h")u' 

^lOf + t4[306u'2 + 90{h^ + h'^)] + t^[288u'^ + GOu'^h^ + h'^) + 30(/i^ + 6h^h'^ + /i'^)] 



1440 

2880 V 

+ 17{h* + h'*) + Q2h^h'^^ + 0(/3'), (36) 

This expansion up to (3"^ is equivalent to the one obtained by Charret et al. Moreover we 
could get two more terms higher than theirs in high-temperature expansion of the grand 
potential. From this expansion we get specific heat and magnetic susceptibility per site at 
half-filled and zero field case: 

+0{(3'), (37) 

X = - + ^- ^- -{3t' + u'')u'(3' + —(3t' + 2u'')t'(3' + (^ + ^ + —)(3^ 
A 2 2 4 6^ ^^ 24^ ^^^8 215^^ 

+0{(3'). (38) 

In figure 1 we plot our numerical results of TBA equations of specific heat at U = 4, /i = 
0, A = U/2. They agree very well with high temperature expansion when (3 < 0.1. Dashed 
line is the numerical calculation of TBA equations by Charret et al. Figure 2 is specific heat 
at [/ = 8. Figure 3 is the susceptibility at U = 4, /i = 0, A = [/ /2. 



9 



U=4 




FIG. 1: Specific lieat at f/ = 4, /i = 0, f//2 — A = 0. Points are our results of TBA calculations 
and dashed line is Charret et al's calculation. Full line is high temperature expansion. 

U=8 




FIG. 2: Specific heat at f/ = 8, /i = 0, f//2 - A = 0. 
IV. DISCUSSION AND CONCLUSION 



In this paper, we show that TBA equations and QTM formulations by Jiittner et 
al give completely the same numerical results and coincide also with the high temper- 
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U=4 




/3 



Figure 3 



FIG. 3: Magnetic susceptibility at f/ = 4, /i = 0, U/2 - A = 0. 

ature expansion. Charret et al' s numerical calculation for TBA equations is not cor- 
rect. Probably the convergence condition is too generous. Our numerical calculation 
is done by use of matliematica 4.1. Source file is available from http://www.issp.u- 
tokyo.ac.jp/labs/theory/mtaka/index.html. 

Kawakami's calculation and Jiittner's calculation are based on Kepler's method. TBA 
equations show very slow convergence in Kepler's method. One needs several tens or several 
hundreds of iterations if one take ric = 6. Jiittner's QTM equation is a bit faster but one 
needs several tens of iterations. Charret's method seems to be based on Kepler's method 
as they reported that the TBA equations converged after 2-3 times iterations. This is too 
short and not reliable. We conclude that the discrepancy between TBA equations with 
QTM equations or high-temperature expansion comes from their inappropriate convergence 
conditions. 

We get very fast convergence if we adopt Newton's method. After 5-6 times of itera- 
tions we get sufficient convergence. Moreover in Newton's method the error decreases with 
acceleration. 

We made numerical programs for TBA and QTM equations. In table I, we observe 
the dependence of the numerical results on ric- As ric increases they converge to the high 
temperature expansion and Jiittner's QTM calculation. Both equations give completely the 
same numerical results. Especially for the grand potential the values coincide with 5 figures 
accuracy. We can conclude that both equations are equivalent, although the mathematical 
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equivalence is not yet proved. 

About the S0(4) symmetry reference [2] did not treat explicitly. But in reference [1] the 
symmetry is treated carefully and the same equations are derived. The thermodynamics 
may not be sensitive to the S0(4) symmetry. In conclusion we can use TBA equations for 
thermodynamic quantities of ID Hubbard model. But one needs some numerical technique 
shown in this paper. 

For the XXZ model it is known that the TBA equations and QTM method give the 
completely the same results numerically [15-18, 20]. Recently, it is shown that TBA equa- 
tions can be derived from quantum transfer matrix formulations for this model [13, 19]. An 
intriguing new simple equation, whic has only one unknown function, was also derived both 
from TBA and QTM [14]. In future, we hope to show the equivalence of two formulations 
for the Hubbard model. 
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-g 


e 


entropy 


C 


X 


nc = 6 


14.96245182 


0.801931104 


1.37643829 


0.01943355 


0.05459858 


Uc = 12 


14.96247427 


0.801886862 


1.37643611 


0.01943778 


0.05465946 


Tic = 24 


14.96247646 


0.801882551 


1.37643590 


0.01943818 


0.05467815 


Tic = 48 


14.96247660 


0.801882261 


1.37643588 


0.01943821 


0.05468324 


Tic = 96 


14.96247662 


0.801882219 


1.37643588 


0.01943821 


0.05468587 


(3 exp. 


14.96246886 


0.801890167 


1.37643590 


0.01943825 


0.05468635 


JKS 


14.96246887 


0.801890163 


1.37643590 


0.01943818 


0.05468632 



TABLE I: Numerical results of TBA equations for grand potential, energy, entropy, specific heat 
and susceptibility at /3 = .1, f/ = 4, /i = U/2 — A = for various values of ric- We put L = 64 and 
compare with high temperature expansion (36) and Jiittner et al's QTM equations. We find that 
ric = 6 is practically sufficient. 

APPENDIX A: i-EXPANSION FOR TBA EQUATIONS 



First, we expand Cq as power series of U' = 4t/U^ 



eo 



1. 



ln2U 



/-I 



1-3,2C(3), 



-2' '2-4' 3 

Then we have t expansion of constant term of — 



{A - eo)/3 = /3A + ^(/^T)^ - ^,m' + 0{fitf. 
We put A = U'x. (Tq^Po^Hq and Si are written as follows: 



^To(A) 



— sech— x 
U 2 



l + ^(-l + 2tanh2^x) + 0(t4) 
(7^ 2 



1 2tln2 , 
Po A; = — + —FT cos k + O , 

Ko{k) = -2tcosk h 0{f ) 

3i(A) = — ^sech— + 



Integral dAs{A — sin A;)... becomes 

2TTt 



dxsix] 



(Al) 



(A2) 



(A3) 



IH sm A; tanh — x + I smA;) ( h tanh —x)-\-0(t 

U 2 ^ U ' ^ 2 2 ' ^ 
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where 



TBA equations are 



s(x) = —sech—x. 
^ ' 4 2 



1 + r]i{x) 



U 4t^ln2 /"^ 

ln({k) = (3{ 2tcosA; — )+ / dxs{x) 

2 U J -oo 

2TTt . TT ,27rt . 2/ 1 T 2 \ 

X IH sm A; tann — x + smA; h tann —x] 

I U 2 ^ U ' ^ 2 2 ' 

lnr]i[x) = — — s[x) — — / s[x)[l-\ — — - sm A; tann — x) m(l + C(A;)) cos A;aA; 

U U J _^ U 2 

+s * ln(l + r]2{x)), 

At r _ 2TTt TT 

\nri'-^{x) = —— I 's{x){l -\ — — - sin A; tanh — x) ln(l + C(A;)) cos A;c/A; + s * ln(l + ?72(^))7 

U J-TT U 2 

\nT]^{x) = s*ln[(l + + 

\ni^{x) = s * ln[(l + + rj'^+,{xj)]- (A4) 

—g(3 is given by 

-gf3 = (3A+ + r + ^ cos k) ln(l + C{k))dk 

pU 2tt ttU 

+ / + -— (-1 + 2tanh-x))ln(l + ?7i(x))c/x. (A5) 

J -oo ^ ^ 

We expand fuctions by power series of 

ln(l + ri,{x)) = In + mf^'\x) + {f3tYff\x) + 



ln(l + rj'^ixj) = ln-^ + {fit)ff'\x) + {fitf ff'\x) 



+ ..., 



3 



ln(l + C{k)) = In ^ + {j3t)z^'Hk) + ^ ____ ^^g^ 

^0 — 1 



We get 



lnry,(x) = In + + (x) + -(a, - oj)if^'\x)y] + 

Oij — i ■' 2 J J 

Inry^x) = + {(3ty,fi'\x) + {(3tY[a'jf\x) + ^ia'^ - af){ff'\x)Y] + 

InCik) = In^- + mzoZ^'Hk) + ifJtY[zoZ^'Hk) + ^(^o - ^o')(^^'H^))'] + (A7) 
zq- 1 2 



where 



[j + 1]^ , _ [j + If _ , ^ ,^/2 [2] 
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and [j] and [j]' are g-integers defined by 

_ s'mh J (3h _ smhj(3{U/2 - A) 

sinh/3/i' ~ smh(3{U/2- A) ' ^ ' 

We have 

jf{x) = f'l'\x) = 0, z^'\k) = -2z-' cos k. (10) 
Equations for f^:^\x), f'^^\x), z^'^^k) are 



(x'ifi^\x) = ^z-'s{x) + s* f!}^\x), 
ajf\x) = s.if^%ix) + f^%ix)), 

a;./f)(x) = 3*(/;(_1(x) + /;g(x)), 

ZoZ^^\k) = 2(1 - ^o"') cos^ k - 

pu 

+ dxs{x){fi'\x)-fi'\x)). (11) 



Same Icind of equations liave appeared in liigli temperature expansion of XXX model. See 
p. 124-126 of [1] or [15]. Solutions of these equations are 



-^'^'W = ^(cosU-^(l + |e-«-)), 

- / X _ J 



Substituting (10) and (12) into (A6) and (A5) we get 

-g(3 = (3A- ln(l - ^o"^) + ln(2 cosh ph) 

+ m'C-^ + i ly'^k) + ^-^ cos kzi'\k)]dk 

+ / 7s{x)[f'\x) + In ^ (-1 + 2tanh^ -x)]t/x| 



(12) 



OO 



= Ine + ^[cosh(/3/.)e^^(l + e^(^^-^)) + -|-e^^^(l - e"^^)], (13) 

where is defined in (33). This coincides with (32) up to the second order. Thus we have 
proven that TBA equations give correct t expansion of —g(3 up to [tfSy. 
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